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Abstract 

Cortical networks, in-vitro as well as in-vivo, can spontaneously generate a 
variety of collective dynamical events such as network spikes, UP and DOWN 
states, global oscillations, and avalanches. Though each of them have been 
variously recognized in previous works as expressions of the excitability of 
the cortical tissue and the associated nonlinear dynamics, a unified picture of 
their determinant factors (dynamical and architectural) is desirable and not 
yet available. Progress has also been partially hindered by the use of a variety 
of statistical measures to define the network events of interest. We propose 
here a common probabilistic definition of network events that, applied to 
the firing activity of cultured neural networks, highlights the co-occurrence 
of network spikes, power-law distributed avalanches, and exponentially dis¬ 
tributed ‘quasi-orbits’, which offer a third type of collective behavior. A 
rate model, including synaptic excitation and inhibition with no imposed 
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topology, synaptic short-term depression, and finite-size noise, acconnts for 
all these different, coexisting phenomena. We find that their emergence is 
largely regulated by the proximity to an oscillatory instability of the dynamics, 
where the non-linear excitable behavior leads to a self-amplification of activity 
fluctuations over a wide range of scales in space and time. In this sense, 
the cultured network dynamics is compatible with an excitation-inhibition 
balance corresponding to a slightly sub-critical regime. Finally, we propose 
and test a method to infer the characteristic time of the fatigue process, 
from the observed time course of the network’s firing rate. Unlike the model, 
possessing a single fatigue mechanism, the cultured network appears to show 
multiple time scales, signalling the possible coexistence of different fatigue 
mechanisms. 


Introduction 


The spontaneous activity of excitable neuronal networks exhibits a spectrum 
of dynamic regimes ranging from quasi-linear, small fluctuations close to 
stationary activity, to dramatic events such as abrupt and transient syn¬ 
chronization. Understanding the underpinnings of such dynamic versatility 
is important, as different spontaneous modes may imply in general differ¬ 
ent state-dependent response properties to incoming stimuli and different 
information processing routes. 

Cultured neuronal networks offer a controllable experimental setting to 
open a window into the network excitability and its dynamics, and have been 
used intensively for the purpose. 

Neuronal cultures in early development phases naturally show alternating 
quasi-quiescent states and ‘network spikes’ (NS) of brief outbreaks of network 
activity 

In addition, recent observations in-vitro (and later even in-vivo) revealed a 
rich structure of network events (‘avalanches’) that attracted much attention 
because their spatial and temporal structure exhibited features (power-law 
distributions) reminiscent of what is observed in a ‘critical state’ of a physical 
system (see e.g. [^|^, and j9, 10 and references therein). Generically, an 
avalanche is a cascade of neural activities clustered in time; while there persist 
ongoing debate on the relation between observed avalanches and whatever 


‘criticality’ may mean for brain dynamics 11 , understanding their dynamical 
origin remains on the agenda. 
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Quasi-synchronous NS, avalanches and small activity flnctnations are 
frequently coexisting elements of the network dynamics. Besides, as we will 
describe in the following, in certain conditions one can recognize network 
events which are clearly distinct from the mentioned network events, which 
we name here as ‘quasi-orbits’. 

The abundant modeling literature on the above dynamical phenomena has 
been frequently focused on specific aspects of one of them 12,13 ; on the other 
hand, getting a unified picture is made often difficult by different assumptions 
on the network’s structure and constitutive elements and, importantly, by 
different methods used to detect the dynamic events of interest. 

In the present paper we define a common probabilistic criterion to detect 
various coexisting dynamic events (NS, avalanches and quasi-orbits) and 
adopt it to analyze the spontaneous activity recorded from both cultured 
networks, and a computational rate model. 

Most theoretical models accounting for NS are based on an interplay 
between network self-excitation on one side, and on the other side some 
fatigue mechanism provoking the extinction of the network spike 


12,13 . For 


such a mechanism two main options, up to details, have been considered: 
neural ‘spike-frequency adaptation’ |^[I4| and synaptic ‘short-term depression’ 
(STD) [^[5,15-18 . Despite their differences, both mechanisms share the basic 


property of generating an activity-dependent self-inhibition in response to 
the upsurge of activity upon the generation of a NS, the more vigorously, the 
stronger the NS (i.e. the higher the average firing rate). In this paper, we 
will mainly focus on STD, stressing the similarities of the two mechanisms, 
yet not denying their possibly different dynamic implications. 

While STD acts as an activity-dependent self-inhibition, the self-excitability 
of the network depends on the balance between synaptic excitation and inhi¬ 
bition; investigating how such balance, experimentally modifiable through 
pharmacology, influences the dynamics of spontaneous NSs is interesting and 
relevant as a step towards the identification of the ‘excitability working point’ 
in the experimental preparation. 

To study the factors governing the co-occurrence of different network 
events and their properties we adopt a rate model for the dynamics of the 
global network activity, that takes into accounts finite-size fluctuations and 
the synaptic interplay between one excitatory and one inhibitory population, 
with excitatory synapses being subject to STD. 

On purpose we implicitly exclude any spatial topology in the model, which 
is meant to describe the dynamics of a randomly connected, sparse network. 
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since we intend to expose the exqnisite implications of the balance between 
synaptic excitation and inhibition, and the activity-dependent self-inhibition 
due to STD. In doing this, we purposely leave out not only the known relevance 
of a topological organization |9, 19,20 , but also the role of cliques of neurons 


which have been proposed to play a pivotal role in the the generation of NS 
as functional hubs 


21 , as well as the putative role of ‘leader neurons’. 


We perform a systematic numerical and analytical study of NSs for varying 
excitation/inhibition balance. The distance from an oscillatory instability of 
the mean-field dynamics (in terms of the dominant eigenvalue of the linearized 
dynamics) largely appears to be the sole variable governing the statistics of the 
inter-NS intervals, ranging from a very sparse, irregular bursting (coefficient 
of variation CV ~ 1) to a sustained, periodic one (CV ~ 0). The intermediate, 
weakly synchronized regime (CV ~ 0.5), in which the experimental cultures 
are often observed to operate, is found in a neighborhood of the instability 
that shrinks as the endogenous fluctuations in the network activity become 
smaller with increasing network size. 

Moreover, the model robustly shows the co-presence of avalanches with 
NS and quasi-orbits. The avalanche sizes are distributed according to a 
power-law over a wide region of the excitation-inhibition plane, although the 
crossing of the instability line is signaled by a bump in the large-size tail of 
the distribution; we compare such distributions and their modulation (as well 
as the distributions of NS) across the instability line with the experimental 
results from cortical neuronal cultures; again the results appear to confirm 
that neuronal cultures operate in close proximity of an instability line. 

Taking advantage of the fact that the sizes of both NS and quasi-orbits 
are found to be significantly correlated with the dynamic variable associated 
with STD (available synaptic resources) just before the onset of the event, we 
developed a simple optimization method to infer, from the recorded activity, 
the characteristic time-scales of the putative fatigue mechanism at work. 
We first tested the method on the model, and then applied it to in-vitro 
recordings; we could identify in several cases one or two long time-scales, 
ranging from few hundreds milliseconds to few seconds. 

Weak or no correlations were found instead between avalanche sizes and 
the STD dynamics; this suggests that avalanches originate from synaptic 
interaction which amplifies a wide spectrum of small fluctuations, and are 
mostly ineffective in eliciting a strong self-inhibition. 
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Models and Analysis 

Experimental data 

As originally described in [^, cortical nenrons were obtained from newborn 
rats within 24 honrs after birth, following standard procednres. Briefly, the 
nenrons were plated directly onto a snbstrate-integrated multielectrode array 
(MEA). The cells were bathed in MEM supplemented with heat-inactivated 
horse serum (5%), glutamine (0.5 mM), glucose (20 mM), and gentamycin 
(10 /ig/ml) and were maintained in an atmosphere of 37° C, 5% C02/95% 
air in a tissue culture incubator as well as during the recording phases. The 
data analyzed here was collected during the third week after plating, thus 
allowing functional and structural maturation of the neurons. MEAs of 60 
Ti/Au/TiN electrodes, 30 /im in diameter, and spaced 200 /im from each other 
(Multi Channel Systems, Reutlingen, Germany) were used. The insulation 
layer (silicon nitride) was pretreated with poly-D-lysine. All experiments were 
conducted under a slow perfusion system with perfusion rates of ~100 pl/h. A 
commercial 60-channel amplifier (B-MEA-1060; Multi Channel Systems) with 
frequency limits of 1-5000 Hz and a gain of 1024 x was used. The B-MEA- 
1060 was connected to MCPPlus variable gain filter amplifiers (Alpha Omega, 
Nazareth, Israel) for additional amplification. Data was digitized using two 
parallel 5200a/526 analog-to-digital boards (Microstar Laboratories, Bellevue, 
WA). Each channel was sampled at a frequency of 24000 Hz and prepared 
for analysis using the AlphaMap interface (Alpha Omega). Thresholds (8x 
root mean square units; typically in the range of 10-20 /iH) were defined 
separately for each of the recording channels before the beginning of the 
experiment. The electrophysiological data is freely accessible for download at 
marom.net.technion.ac.il/neural-activity-data/. 


Network rate dynamics 


A set of Wilson-Cowan-like equations 22 for the spike-rate of the excitatory 


(z/^;) and the inhibitory (z//) neuronal populations lies at the core of our 
dynamic mean-field model: 


— —{^E — ^{Ie)) ( 1 ) 

TiUi = -(z//- <!)(//)), 
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where te and tj represent two characteristic times (of the order of few to few 
tens of ms), and $ is the gain function of the input currents, Ie and //, that 
in turn depend on ue, and the synaptic efficacies. We chose <J> to be the 
transfer function of the leaky integrate-and-fire neuron under the assumptions 
of Gaussian, uncorrelated input of mean /i and infinitesimal variance 23 


$[/r, cr^] = y/wTv 






f yreset _Yrest 


exp (s^) [erf(s) + l] ds + Trefract 


-1 


\A 


ry 


( 2 ) 

where Ty is the membrane time constant, Trefract is a refractory period, and 
yrest^ yreset^ ythresh respectively the rest, the post-firing reset, and the 

firing-threshold membrane potential of the neuron (we assume the membrane 
resistance R = 1). 

The model incorporates the non-instantaneous nature of synaptic trans¬ 
mission in its simplest form, by letting the vs being low-pass filtered by a 
single synaptic time-scale f: 


T V = {v — v). (3) 

One can regard the variables vs as the instantaneous firing rates as seen by 
post-synaptic neurons, after synaptic filtering. The form of Eq. and our 
choice of f values (see Table implicitly neglects slow NMDA contributions 
and is restricted to AMPA and GABA synaptic currents. Thus, the input 
currents Ie and // in Eq. [^will be functions of the rates vs through these 
filtered rates; with reference to Eq. the model assumes the following form 
for the mean and the variance of the current Ie (the expressions for 1/ are 
similarly defined): 

= CHe Ve Wexc JeE Te + 

C ni Vi Winh JeI + I'ext Jext (4) 

Ct| = CnEVE^l^e{JEE + ^%E)^E + 

cni Vi wf„h {Jei + {Jjxt + > 

where the Ue and ni are the number of neurons in the excitatory and inhibitory 
population respectively; c is the probability of two neurons being synaptically 
connected; Jee {Jei) is the average synaptic efficacy from an excitatory 
(inhibitory) pre-synaptic neuron to an excitatory one, cxj is the variance 
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of the J-distribution; Wgxc and Wjnh are dimensionless parameters that we 
will nse in the following to independently rescale excitatory and inhibitory 
synapses respectively. Finally, an external current is assumed in the form of 
a Poisson train of spikes of rate Vext driving the neurons in the network with 
average synaptic efficacy Jext- In Eq- (0 < < 1) is the fraction 

of synaptic resources available at time t for the response of an excitatory 
synapse to a pre-synaptic spike; the evolution of evolves according to the 
following dynamics, which implements the effects of short-term depression 
(STD) (M[^ 


into the network dynamics: 


"Tstd f'E — (1 — 'Te)— mstd rE t'std 


(5) 


where 0 < mstd < 1 represents the (constant) fraction of the available synaptic 
resources consumed by an excitatory postsynaptic potential, and tstd is the 
recovery time for the synaptic resources. 

Finally, for a network of n neurons, we introduce finite-size noise by 
assuming that the signal the synapses integrate in Eq. is a random process 
Un of mean u; in a time-bin dt, we expect the number of action potentials 
fired to be a Poisson variable of mean nh>{t) dt; Eq. |^will thus become: 

TV = {Vn-v) 

_ Poisson [n i/dt] 


Putting all together, the noisy dynamic mean-field model is described by 


(stochastic) differential equations: 

TeVe 

= <h(/i£;, cr|) - i/ij 

TjVj 

= (rfj - Vi 

TeVe 

= VnE - Ve 

fiVl 

= Vnj - Vl 

'Tstd te 

= — Te) — UstdTsTBTeVe 


(7) 


complemented by Eqs. ii and The values of all the fixed network 
parameters are shown in Table Since we will compare the dynamics of 
networks of different sizes, we scale the connectivity with network size in 
order to keep invariant the mean field equations: we hold the number of 
synaptic connection per neuron constant by rescaling, with reference to Eq. 
the probability of connection c so that cue and cri/ are kept constant to the 
reference values that can be deduced from Table [H 
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Table 1. Network parameters. 


Parameter 

Value 

He / nj 

160 / 40 

'^V / ^refract 

20 / 2 ms 

'^rest j firing 

-70 / -55 mV 

c 

0.25 

JeE i ^Jee 

0.809 ± 0.202 mV 

Jj E i ^ Ji E 

1.23 ± 0.307 mV 

J El ^ ^ Jei 

-0.340 ± 0.0850 mV 

Jii dz ajjj 

-0.358 ± 0.0894 mV 

Jext i 

0.416 ±0.104 mV 

^ext 

1.25 kHz 

Te / Ti 

10 / 2 ms 

Tr 

800 ms 

USTD 

0.2 

Te = Ti 

20 ms 


Spike-frequency adaptation (SFA) (not present in simulations unless where 
explicitly stated) is introduced by subtracting a term to the instantaneous 
mean value of the Ie current: 


flE ^ — 9SFA CE{t) 


( 8 ) 


proportional to the instantaneous value of the variable cg, that simply inte¬ 
grates VriE- 

t-sfa = -Ce + (9) 

with a characteristic time tsfa- This additional term aims to model an after¬ 
hyperpolarization, Ca^"*"-dependent K+ current 26,2m. In this sense, ce can 


be interpreted as the cytoplasmic calcium concentration [Ca^"'']), whose effects 
on the network dynamics are controlled by the value of the “conductance” 
5'sfa- 

Simulations are performed by integrating the stochastic dynamics with a 
fixed time step dt = 0.25 ms. 

In the following, by “spike count” we will mean the quantity u{t)n dt. 





















Network events detection 


For the detection of network events (NSs, quasi-orbits, and avalanches) we 
developed a unified approach based on Hidden Markov Models (HMM) |^ . 
Despite HMM have been widely used for temporal pattern recognition in many 
different fields, to our knowledge few attempts have been made to use them in 
the context of interest here 29,30 . For the purpose of the present description, 


we just remind that a HMM is a stochastic system that evolves according 
to Markov transitions between “hidden”, i.e. unobservable, states; at each 
step of the dynamics the visible output depends probabilistically on the 
current hidden state. Such models can be naturally adapted to the detection 
of network events, the observations being the number of detected spikes 
per time bin, and the underlying hidden states, between which the system 
spontaneously alternates, being associated with high or low network activity 
(‘network event - no network event’). A standard optimization procedure 
adapts then the HMM to the recorded activity sample by determining the 
most probable sequence of hidden states given the observations. 

The two-step method we propose is based on HMM, has no user-defined 
parameters, and automatically adapts to different conditions. 

In the first step, the algorithm finds the parameters of the two-state HMM 
(one low-activity state, representing the quasi-quiescent periods, and one 
high-activity state, associated with network events) that best accounts for a 
given sequence of spike counts - the visible states in the HMM; such fitting 
is performed through the Baum-Welch algorithm 28 . In the second step. 


the most probable sequence for the two alternating hidden levels, given the 
sequence of spike counts and the fitted parameters, is found through the 
Viterbi algorithm. Network events are identified as the periods of dominance 
of the high activity hidden state. 

In order to retain only the most significant events a minimum event 
duration is imposed; such threshold is self-consistently determined as follows. 
The Viterbi algorithm is also applied to a “surrogate” time-series obtained by 
randomly shuffling the original one, thereby generating a set of “surrogate” 
events. The purpose is to determine the desired minimum event duration 
from the high duration tail of surrogate events (which, by construction, 
come from a time-series with no real temporal structure). Since the high 
duration distribution tail is found to be roughly exponential, we fit such tail 
by considering only the surrogate events of duration larger than the 75th 
percentile. Then, from the fitted exponential, we compute the duration value 
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such that the probability of durations greater than this value is P(surrogate) = 
10“^. In other words, we set the threshold on minimum duration of detected 
events to the duration of exceptionally long (P < 10“^) surrogate events. 

As already remarked, we used essentially the same algorithm for detecting 
NS/quasi-orbits and avalanches. The only significant difference is that, in 
the case of avalanches, the emission probability of the low-activity hidden 
state is kept fixed during the Baum-Welch algorithm to pin) ~ hnO (hij 
is the Kronecker delta; p{n) is the probability of emitting n spikes in a 
time-bin). Thus the lower state is constrained to a negligible probability of 
outputting non-zero spike-counts, conforming to the intuition that in between 
avalanches the network is (almost) completely silent. More precisely, we 
set p{l) = 10“®(n), where (n) is the average number of spikes that the 
network emits during a time-bin dt. After the modified Baum-Welch first 
step, avalanches are determined, as above, by applying the Viterbi algorithm; 
no threshold is applied in this case, neither to the avalanche duration nor to 
its size. 

The proposed procedures introduce three arbitrary parameters: the time 
bin dt, the probability P(surrogate) for network spikes and quasi-orbits, and 
the probability p(l). To test the robustness of the algorithms, we varied 
these parameters over ample ranges: dt between 0.25 and 8 ms; P(surrogate) 
between 10“^ and 10“"^; p{l) between 10“® and 10“^. We found that avalanche 
size distributions are virtually unaffected under variations of p(l), and only 
mildly affected for the largest dt explored; higher values of P(surrogate) 
lead, as expected, to detect a larger number of small quasi-orbits, yet these 
additional events do not alter the overall shape of the size distribution 
predicted by the theory (see next section); on the other hand, a large number 
of very small quasi-orbits does have a detrimental effect on the correlation 
results reported in Section “Inferring the time-scales”. 

Simulations and data analysis have been performed using custom-written 
mixed C++/ MATLAB (version R2013a, Mathworks, Natick, MA) functions 
and scripts. 

Size distribution for quasi-orbits and network spikes 

The non-linear rate model described above can show a wide repertoire of 
dynamical patterns, as for example multiple stable fixed points and large, 
quasi-periodic oscillations. As we will show, for sufficiently excitable networks, 
a stable state of asynchronous activity (fixed point) is destabilized, in favor 
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of stable global oscillations. Finite size noise probes differently network’s ex¬ 
citability at different distances from such instability. Before global oscillations 
become stable (in the infinite network limit), the network’s highly non-linear 
reaction to its own fluctuations can ignite large, relatively stereotyped, “net¬ 
work spikes”. Also, in the proximity of the oscillatory (Hopf) instability, noise 
can promote “quasi orbits”, he., transient departures from the fixed point 
which develop on time-scales dictated by the upcoming oscillatory instability, 
of which they are precursors. Under a linear approximation, the probability 
distribution of the amplitude I of these quasi-orbits can be explicitly derived 
as explained in the following. 

Consider a generic planar linear dynamics with noise: 

z = At. + a (10) 

where ^ is 2X 2 real matrix, and ^ = (^(t), 0) is a white noise with = 

S{t—t'). We here assume that the system is close to a Hopf bifurcation; in other 
words that the matrix A has complex-conjugated eigenvalues \± = 3f?A -f i AA, 
with 3f?A < 0 and |3fJA| AA. 

By means of a linear transformation, the system can be rewritten as: 

X = 3f?A X — AA y + ax ^ 

y = ^Xx + ^Xy + ay^, 

(11) 

with ax and ay constants determined by the coordinate transformation. Mak¬ 
ing use of Itb’s lemma to write: 

x"^ = 2'SiXx‘^ — 2‘^Xxy + al + 2x ax ^ 
y^ = 2^Xxy + 2^Xy^ + al + 2yay^, 

and summing the previous two equations, we find for the square radius 
1“^ = x"^ + y'^ the dynamics: 

P = 2 3f?A P -|- -|- 2 (^ax x + ayy) (12) 


with = al + ay. 

As long as AA 3> |3^A|, it is physically sound to make the approximation: 

{x{t), y{t)) = /(O) ( cos(AA t + (p), sin(AA t + 0)), (13) 
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for 0 < t < T = 2 tt/Q'A and then to average the variance of the noise over 
such period to get: 




T 


= m) 


2 + 


cos(3A t + (j)) + ay sin(^rA t + 0) 

/(0)2a'2 


dt = 


in order to rewrite Eq. (12) as: 


P = 2mi^ + a'^ + V2la'^. 


(14) 


Such stochastic differential equation is associated with the Fokker-Planck 
equation: 


dtp{l‘^,t) = —di2 [2 3fJA/^ + a'‘^]p{P,t) + 
+ Cr'^ d'^2 p{f, t) = L;2 p{P, t) 


( 16 ) 


that admits an exponential distribution as stationary solution: 


Pssin = 


2, 2|3fJA| . 2|3fJA|/2 


a‘ 


/2 


exp 


a‘ 


/2 


). 


that is, a Rayleigh distribution for /: 


4|gfJA|, , 2|gfJA|/% 

Pss{l) = —;j^lexp{ - 


a'- 


( 16 ) 


(17) 


On the other hand, we found a correlation between I (the maximal depar¬ 
ture from the low-activity fixed point) and the duration of the quasi-orbit. 
Therefore the size of the quasi-orbit (the ‘area’ below the firing rate time 
profile during the excursion from the fixed point) is expected to scale as P, 
so that it should be exponentially distributed. 

For network spikes we do not have a theoretical argument to predict the 
shape of the size distribution, however empirically a (left-truncated) Gaussian 
distribution proved to be roughly adequate. Since we expect that quasi-orbits 
and NS contribute with different weights for varying excitatory/inhibitory 
balance, we adopted the following form for the overall distribution of network 
event size to fit experimental data: 


p{x) 


^0 


exp ( 


(X - Xq) 

To 


) + 


1 -Po 

\/^ ai 


exp ( 


(x — 

2af ^ 


(18) 
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The parameters of the two distributions and their relative weight 0 < po < 1 
are estimated by minimizing the log-likelihood on the data. A threshold 
for the event size is determined as the value having equal probability of 
being generated by either the exponential or the normal distribution. In the 
following, NSs are defined as events having size larger than this threshold. 
In those cases in which a threshold smaller than the peak of the normal 
distribution could not be determined, no threshold was set. 


Results 


In the following, we will study a stochastic firing-rate model and make 
extensive comparison of its dynamical behavior with the activity of ex-vivo 
networks of cortical neurons recorded through a 60-channel multielectrode 
array. 

The first question we want to answer is how the excitation-inhibition 
balance affects network dynamics. Starting from the statistics of network 
spikes (NS) we show that it is well described by a single variable measuring 
the distance from an oscillatory instability of the dynamics. We then study 
in the model the effects of finite-size fluctuations on the statistics of NS. 

Then, taking advantage of new detection algorithm we introduce (see 
Models and Analysis), we recognize the presence of a spectrnm of network 
events, including three families: NS, “quasi-orbits”, and avalanches. The 


predicted size distribution of quasi-orbits , exponential component in Eq. 18 


is confirmed by simulations and recovered in experimental data analysis. We 
investigate how the different network events characterize in various propor¬ 
tions the network dynamics depending on the excitatory-inhibitory balance; 
experimental data offer an interesting match with model findings, compatible 
with ex-vivo network being typically slightly below the oscillatory instability. 

Finally we introduce a simple procedure to infer the time-scales of putative 
slow self-inhibitory mechanisms underlying the occurrence of network events. 
The inference is obtained based on knowledge of the firing activity alone; this 
makes the method interesting for analysis of experimental data, as we show 
through exemplary results. 

The stochastic firing-rate model consists of two populations of nenrons, 
one excitatory and one inhibitory, interacting through effective synaptic 
couplings; excitatory synaptic couplings follow a dynamics mimicking short¬ 
term depression (described after the Tsodyks-Markram model, (^). We 
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adopted the transfer function of the leaky integrate-and-fire neuron subject 
to white-noise current with drift 23 as the single population input-output 
function; moreover the activity of each population is made stochastic by 
adding realistic finite-size noise. Working with a noisy mean field model 
allows in principle to easily sweep through widely different network sizes and, 
more importantly, allows us to perform the stability analysis. 

To start the exploration that follows, we chose a reference working point 
where the model’s dynamics has a low-rate fixed point (2 — 4 Hz) just on 
the brink of an oscillatory instability or, in other words, where the dominant 
eigenvalue A of the dynamics, linearized around the fixed point, is complex 
with null real part. The model network (Fig.[^ panel A) shows in proximity of 
this point a dynamical behavior qualitatively similar, in terms of population 
spikes, to what is observed in ex-vivo neuronal networks (Fig. panel B). 



t (s) 

Figure 1. Time course of the network firing rate. Panel A: noisy 
mean-field simulations; panel B: ex-vivo data. Random large excursions of 
the firing rate (network spikes and quasi-orbits) are clearly visible in both 
cases. 
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Excitation-inhibition balance and network spike statis¬ 
tics 

As the relative balance of excitation and inhibition is expected to be a 
major determinant of NS statistics we investigated first, for spontaneous NSs, 
how the inter-NS intervals (INSI) and their regularity (as measured by the 
coefficient of variation, CVinsi) depend on such balance. In Fig. [^we report 
the average INSI (left panel) and CVinsi (right panel) in the plane (wexc, 
Winh) of the excitatory and inhibitory synaptic efficacies {Jee —t weJee, 
JiE —)■ We Jie, Jei —t wi Jei, Jii —t wi Jji, see Eq. i- Starting from the 
center of this plane (wexc = 1, Winh = 1) and moving along the horizontal axis, 
all the excitatory synapses of the network are multiplied by a factor Wexc^ 
moving right, the total excitation of the network increases (wexc > 1), toward 
left it decreases (wexc < !)• Along the vertical line, instead, all the inhibitory 
synapses are damped (moving downward, Winh < 1) or strengthened (going 
upward, winh > !)• 

It is clearly seen that both (INSI) and CVinsi are approximately dis¬ 
tributed in the plane along almost straight lines of equal values: for a chosen 
(INSI) or CVinsi one can trade more excitation for less inhibition keeping 
the value constant, suggesting that, at this level of approximation, a measure 
of net synaptic excitation governs the NS statistics. Besides, not surprisingly, 
for high net excitation NSs are more frequent (~ 1 Hz) and quasi-periodic 
(low CVinsi), due to the fact that the STD recovery time determines quasi- 
deterministically when the network is again in the condition of generating a 
new NS. Weak excitability, on the other hand, leads to rare NSs, approaching 
a Poisson statistics (CVinsi — 1), since excitability is so low that fluctuations 
are essential for recruiting enough activation to elicit a NS, with STD playing 
little or no role at the ignition time; below an “excitation threshold”, NSs 
disappear. 

The solid lines in Fig. are derived from the linearization of the 5- 
dimensional dynamical system (see Eq. [^, and are curves of iso-jPA, where A 
is the dominant eigenvalue of the Jacobian: JJA = 0 Hz (white line, signaling 
a Hopf bifurcation in the corresponding deterministic system), JJA = 3.5 Hz 
(red line), and JJA = —3.5 Hz (black line). Values of CV found in typical 
cultured networks are close to model results near the bifurcation line JJA = 0 
Hz. We observe, furthermore, that such lines roughly follow iso-(INSI) and 
iso-CViNsi curves, suggesting that a quasi one-dimensional representation 
might be extracted. 
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Figure 2. Inter-network-spike interval (INSI) statistics in the noisy 
mean-field model, for varying levels of excitation (wexc) ^nd inhibi¬ 
tion (winh)- Panel A: (INSI) (the scale is in seconds); panel B: coefficient of 
variation of INSI (CVinsi)- For high net excitation (bottom-right quadrant) 
short-term depression plays a determinant role in generating frequent and 
regular (low CVinsi) NSs; for weak excitability (upper-left quadrant) random 
fluctuations are essential for the generation of rare, quasi-Poissonian NSs 
(CVinsi — !)• The solid lines are isolines of the real part IPA of the dominant 
eigenvalue of the mean-field dynamics’ Jacobian; white line: 3fJA = 0 Hz; red 
line; IRA = 3.5 Hz; black line; 3RA = —3.5 Hz. Note how such lines roughly 
follow isolines of (INSI) and CVinsi- 


We show in Fig. (INSI) (panel A) and CVinsi (panel B) against 3RA 
for the same networks (circles) of Fig. and for a set of larger networks 
{N = 8000 neurons, squares) that are otherwise identical to the first ones, 
pointwise in the excitation-inhibition plane (the average number of synaptic 
connections per neuron for the larger networks is kept constant to the value 
used in the original, smaller ones, as explained in Models and Analysis) The 
difference in size amounts, for the new, larger networks, to weaker endogenous 
noise entering the stochastic dynamics of the populations’ firing rates (see 
Eq. second line). The points are seen to approximately collapse onto 
lines for both sets of networks, thus confirming 3RA as the relevant control 
quantity for (INSI) and CVinsi- It is seen that, for the smaller networks, 
(INSI) and CVinsi depend smoothly on 3RA, due to finite-size effects smearing 
the bifurcation. Also note the branch of points (filled circles) for which 
A"A = 0 and then no oscillatory component is present, corresponding to points 
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in the extreme top-left region of the planes in Fig. For the set of larger 
networks, the dependence of (INSI) and CVinsi on the 3fJA is much sharper, as 
expected given the much smaller finite-size effects; this shrinks the available 
region, around the instability line, allowing for intermediate, more biologically 
plausible values of CVinsi- 




Figure 3. Stability analysis of the linearized dynamics captures 
most of the variability in the inter-network-spike interval (INSI) 
statistics. (INSI) (panel A) and CVinsi (panel B) vs the real part 3f?A of 
the dominant eigenvalue of the Jacobian of the linearized dynamics, for two 
networks that are pointwise identical in the excitation-inhibition plane, except 
for their size (circles: 200 neurons, as in Fig. squares: 8000 neurons). The 
data points almost collapse on 1-D curves when plotted as functions of JJA, 
leading effectively to a “quasi one-dimensional” representation of the INSI 
statistics in the (wexc, ■Winh)-plane. The region in which the INSIs are neither 
regular (CVinsi 0) nor completely random (CVinsi — 1), as typically 
observed in experimental data, shrinks for larger networks. The filled circles 
mark a null imaginary part Q'A. 


We remark that NSs are highly non-linear and relatively stereotyped 
events, typical of an excitable non-linear system. The good predictive power 
of the linear analysis for the statistics of INSI signals that relatively small 
fluctuations around the system’s fixed point, described well by a linear analysis, 
can ignite a NS. 
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A spectrum of network events 


Our mean-field, finite-size network is a non-linear excitable system which, 
to the left of the Hopf bifurcation line, and close to it, can express different 
types of excnrsions from the otherwise stable fixed point. Large (almost 
stereotyped for high excitation) NSs are exquisite manifestations of the non¬ 
linear excitable nature of the system, ignited by noise; the distribution of NS 
size (nnmber of spikes generated dnring the event) is relatively narrow and 


approximately symmetric (the Gaussian component of Eq. 18). 

Noise can also induce smaller, transient excursions from the fixed point 
which can be adequately described as quasi-orbits in a linear approximation. 
In fact, noise indnces a probability distribntion on the size of snch events, 
which can be computed as explained in Methods and Analysis (the exponential 
part in Eq. 18). Eig. panel A, shows the activity of a simnlated network 
(blue line) alongside with detected network events. We remark that the the 
different event types may not in general be easily distingnished on a single¬ 
event basis, while we argue that they are probabilistically distinguishable. 
From the best fit for the expected size distribution a threshold for the event 
size can be determined to separate events that are {a-posteriori) more probably 
quasi-orbits from the ones that are more probably NSs (for details, see Models 
and Analysis). Following such classification, the green line in Fig. panel A, 
marks the detection of two NSs (first and third event) and two quasi-orbits 
(second and fourth event). 

We also emphasize that the existence of qnasi-orbits is a specific conse¬ 
quence of the fact that in the whole excitation-inhibition plane explored for 
the model, the low-activity fixed point becomes unstable via a Hopf bifur¬ 
cation. It is indeed known that for nonlinear systems in the proximity of a 
Hopf bifurcation, noise promotes precursors of the bifurcation, which appear 
as transient synchronization events (see, e.g., |^). 

As one moves around the excitation-inhibition plane, to the left of the 
bifurcation line, the two types of events contribute differently to the overall 
distribution of network event sizes. Qualitatively, the farther from the bifur¬ 
cation line, the higher the contribntion of the small, “quasi-linear” events. 
This fact can be understood by noting that the average size of such events is 
expected to scale as l/|3fJA|, where 3fJA is the real part of the dominant eigen¬ 


value of the (stable) linearized dynamics (see Models and Analysis, Eq. 16). 


The average size is fnrthermore expected to scale with the amonnt of noise 
affecting the dynamics, thus the contribution of qnasi-linear events is also 
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Figure 4. Algorithms for network events detection. Panel A: total 
network activity from simnlation (blue line) with detected NS/quasi-orbits 
(green line) and avalanches (red line). Four large events (green line) are 
visible; the first and third are statistically classified as network spikes; the 
other smaller two as quasi-orbits. Note how network spikes and quasi-orbits 
are typically included inside a single avalanche. Panel B: a zoom over 0.5 
seconds of activity, with discretization time-step 0.25 ms, illustrates avalanches 
structure (red line). 


expected to vanish for larger networks. 

It has been previously reported that activity dynamics may be different 
from one network to the other, reflecting idiosyncrasies of composition and 
history-dependent processes ( j^). Moreover, the dynamics of a given 
network, as well as its individual neurons, may shift over time (minutes and 
hours) between different modes of activity ( |32|-[34]). We therefore chose 
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to demonstrate the efficacy of our analytical approach on two data sets of 
large-scale random cortical networks. 

In panels A-C of Fig|^ we show the experimental distributions of event sizes 
for two cultured networks: panels A and B are ~40-minute recordings taken 
from a very long recording for the same network; panel C is ~ 1-hour recording 
from a different cultured network. By visual inspection, the distributions 
appear to be consistent with two components contributing with various 
weights, both for different periods of the same network, and for different 
networks. In the light of the above theoretical considerations, one is led to 
generate the hypothesis that the two components contributing to the overall 
distribution were associated with quasi-orbits and network spikes respectively; 
to test this hypothesis, we fitted (solid lines in Fig. the experimental 
distributions with the sum of an exponential and a Gaussian distribution (see 
Models and Analysis, Eq. 18), prepared to interpret a predominance of the 
exponential (Gaussian) component as a lesser (greater) excitability of the 
network. We remark that (see panels A and B) the relative weights of the 
two components appear to change over time for the same network, as if the 
excitability level would dynamically change; more on this at the end of this 
section. 

To substantiate the above interpretation of experimental results, we turned 
to long simulations (about 5.5 hours) of networks in different points in the 
excitation-inhibition plane (Fig. [^, from which we extracted the distribution 
of network events and fitted them with Eq. 18 as for experimental data (see 
panels D-F in Fig. [^. Again, to the eye, the fits appear to be consistent 
with the two components variously contributing to the overall distribution, 
depending on the excitability of the network. 

If, however, the fits are subject to a Kolmogorov-Smirnov test, the test 
fails (p < 0.01) for panels D and F. By inspecting the maximum distance 
between the cumulative distributions for simulation data and the fit, we found 
it at the lowest size bin for panel D, while the “Gaussian” part gives the 
greater mismatch for panel F. As for panel D, while the theoretical argument 
for the quasi-orbits clearly captures the shape of the size distributions, the 
way the test fails in the exponential part is interesting. 

In fact, network events cannot be detected with arbitrarily small size: in 
a way, the detection procedure imposes a soft threshold on the event size, 
below which the exponential distribution is not applicable. 

We can provide a rough estimate of such soft threshold as follows. A 
quasi-orbit duration is, to a first approximation, proportional to 1/AA, which 
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is of the order of few hundreds milliseconds not too far from the bifurcation 
line in the excitation-inhibition plane. Taking, for instance, 150 ms, an 
event will be detected if network activity within this time-span is larger than 
average (typically few spikes per second per neuron; we take 3 for the present 
example): this leads to a soft threshold of abont 100 spikes. This wonld be 
the lower limit of applicability of the exponential part of the distribntion; 
this also explains the trongh observed for very small sizes. 

As for the failure of the Kolmogorov-Smirnov test for the right part of 
the distribution in panel F, it should be remarked that the assumption of 
a Ganssian distribution for the size of network spikes, although generically 
plausible, is not gronnded in a theoretical argnment, and it’s not surprising 
that, on the order of 10^ detected events, even a moderate skewness, as the 
one observed, can make the test fail. 

The fit for experimental data of panels A-C passed the Kolmogorov- 
Smirnov test (p > 0.01). 

As mentioned in the introdnction, avalanches are cascades of nenral 
activities clnstered in time (see Models and Analysis for our operational 
definition; examples of different methods used in the literatnre to detect 
avalanches can be found in [7,35-^). Fig. panel A and panel B, shows an 
example of the strncture of the detected avalanches (red lines) in the model 
network. 

We extracted avalanches from simulated data, as well as from experimental 
data. For simnlations, we choose data corresponding to three points in the 
(Wexc, Winh) plane of Fig. i with constant Wjnh = 1 and increasing Wexc, with 
the rightmost falling exactly over the instability line (white solid line in Fig.[^. 
Three experimental data sets were extracted from different periods of a very 
long recording of spontaneous activity from a neural culture; each data set is 
a 40-minute recording. 

In Fig. [^we show (in log-log scale) the distribution of avalanche sizes for 
the three simulated networks (top row) and the three experimental (bottom 
row) data sets (bine dots); red lines are power-law fits (M) . 

From the panels in the top row we see that the distributions are well 
fitted, over a range of two orders of magnitude, by power-laws with exponents 
ranging from about 1.5 to about 2.2, consistent with the results found in [^. 
Note that in the cited paper the algorithm nsed for avalanche detection 
is quite different from onrs, and the wide range of power-law exponents 
is related to their dependence on the time-window nsed to discretize data. 
In (adopting yet another algorithm for avalanche detection), both the 
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Figure 5. A broad spectrum of synchronous network events: simu¬ 
lations vs ex-vivo data. Panels A-C: experimental distributions of network 
events. Panels A and B; ~40-minute recordings from a very long recording, 
for the same network; panel C: ~ 1-hour recording from another cultured net¬ 
work. Panels D-F: distributions from simulations of networks corresponding 
to the points in Fig. ((Wexc, = (0.82, 0.7), (Wexc, Wi,h) = (0.82, 0.55), 
(wexc, Winh) = (0.88, 0.55)). The three networks of panels D-F have increasing 
levels of subcritical excitability. Note the logarithmic scale on the y-axis. 
The solid lines are fits of the theoretical distribution of event sizes, a sum 
of an exponential (for quasi-orbits) and a Gaussian (for NS) distribution 


(see Models and Analysis, Eq. 18). The vertical lines mark the probabilistic 


threshold separating NS and quasi-orbits. 


shape of the avalanche distribution and the exponent vary depending on using 
pharmacology to manipulate synaptic transmission, over a range compatible 
with our model findings; notably, they find the slope of the power-law to be 
increasing with the excitability of the network, which is consistent with our 
modeling results. 

Panels B and C of Fig. [^clearly show the buildup of ‘bumps’ in the high- 
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Figure 6. Avalanche size distribution: simulations vs ex-vivo data. 

Panels A-C: mean-field simnlations, with fixed inhibition Wjnh = 1- and 
increasing excitation (wgxc = 0.9, 0.94, 1). The distribntions are well fitted 
by power-laws; panel B and C clearly show the bnildnp of ‘bnmps’ in the 
high-size tails, reflecting the increasing contribntion from network spikes 
and quasi-orbits in that region of the distribution. Panels D-F from ex-vivo 
data, different ~ 40-minute segments from one long recording; power-laws 
are again observed, although fitted exponents cover a smaller range; in 
panels E and F, bumps are visible, similar to model findings. The similarity 
between the theoretical and experimental distributions could reflect changes 
of excitatory/inhibitory balance in time in the experimental preparation. 
Since all the three simulations lay on the left of or just on the bifurcation line 
(white line in Fig. [^, the shown results are compatible with the experimental 
network operating in a slightly sub-critical regime. 


size tails, increasing with the self-excitation of the network; this is understood 
as reflecting the predominance of a contribution from NS and possibly quasi¬ 
orbits in that region of the distribution, on top of a persisting wide spectrum 
of avalanches. This feature also is consistent with the experimental findings 
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of 39 , and has been previously shown in a theoretical model l40] for non-leaky 


integrate-and-fire neurons endowed with STD and synaptic facilitation. 

Turning to the plots in the bottom row of Fig. we observe the following 
features: power-laws are again observed over two decades and more; in panels 
E and F, bumps are visible, similar to model findings; power-law exponents 
cover a smaller range just above 2. 

While the sequence of plots in two rows (modeling and experiment) 
clearly shows similar features, we emphasize that experimental data were 
extracted from a unique long recording, with no intervening pharmacological 
manipulations affecting synaptic transmission; on the other hand, it has 
been suggested 41 that a dynamic modulation of the excitatory/inhibitory 
balance can indeed be observed in long recordings; although our model would 
be inherently unable to capture such effects, it is tempting to interpret the 
suggestive similarity between the theoretical and experimental distributions 
in Fig. [^as a manifestation of such changes of excitatory/inhibitory balance 
in time, of which the theoretical distributions would be a ‘static’ analog. 
To rule out the possibility that different behaviors in time could be due 
to intrinsic and global modifications in the experimental preparation, we 
checked (see Fig. ?? - SI Fig) the waveforms of the recorded spikes across all 
MEA electrodes, comparing the earliest and latest used recordings (about 40 
minutes each, separated by about 34 hours). In most cases the waveforms 
for the two recordings are remarkably similar, and when they are not, no 
systematic trend in the differences is observed. 

If our interpretation is correct, the experimental preparation operates 
below, and close, to an oscillatory instability; on the other hand, contrary 
to NS, the appearance of avalanches does not seem to be exquisitely related 
to a Hopf bifurcation, rather they seem to generically reflect the non-linear 
amplification of spontaneous fluctuations around an almost unstable fixed 
point - a related point will be mentioned in the next section. We also remark 
that we obtain power-law distributed avalanches in a (noisy) mean-field rate 
model, by definition lacking any spatial structure; while the latter could well 
determine specific (possibly repeating) patterns of activations (as observed 
in [^), it is here suggested to be not necessary for power-law distributed 
avalanches. 

The avalanche size distribution for the same network as in Fig. panel C, 
is sparser but qualitatively compatible with the distribution in Fig. panel 
F (see Fig. ?? - S2 Fig); in particular, the distribution shows a prominent 
peak for high-size avalanches, consistently with the interpretation, given in 
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connection with Fig. of high excitability. 

We do not provide examples of avalanche and NS-quasi orbits size dis¬ 
tributions in the super-critical region on the right of the Hopf bifurcation 
line in Fig. this is because the phenomenology in that region is relatively 
stereotyped and easy to guess/understand: the high excitability of the network 
generates, moving on the right of the bifurcation line, increasingly stereo¬ 
typed network spikes, which dominate the size distribution of the network 
events (see Fig. ?? - S3 Fig, panel A); even though finite-size fluctuations 
blur the bifurcation line, quasi-orbits are expected to contribute very little 
in the supercritical region; the distribution of avalanche sizes is increasingly 
dominated by the high-size bump associated with network spikes (see Fig. ?? 
- S3 Fig, panel B). 

Inferring the time-scales 

The fatigue mechanism at work (STD in our case) is a key element of the 
transient network events, in its interplay with the excitability of the system. 
While the latter can be manipulated through pharmacology, STD itself (or 
spike frequency adaptation, another neural fatigue mechanism) cannot be 
directly modulated. It is therefore interesting to explore ways to infer relevant 
properties of such fatigue mechanisms from the experimentally accessible 
information, i.e. the firing activity of the network. We focus in the following 
on deriving the effective (activity-dependent) time scale of STD from the 
sampled firing history. 

The starting point is the expectation that the fatigue level just before a 
NS should affect the strength of the subsequent NS. We therefore measured 
the correlation between r (fraction of available synaptic resources) and the 
total number of spikes emitted during the NS (NS size) from simulations. We 
found that the average value of r just before a NS is an effective predictor of 
the NS size, the more so as the excitability of the network grows. 

Based on the r-NS size correlation, we took the above “experimental” 
point of view, that only the firing activity v is directly observable, while r is 
not experimentally accessible. Furthermore, the success of the linear analysis 
for the inter-NS interval statistics (due to the NS being a low-threshold very 
non-linear phenomenon), suggests that without assuming a specific form for 
the dynamics of the fatigue variable /, we may tentatively adopt for it a 
generic linear integrator form, of which we want to infer the characteristic 
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time-scale t*\ 


(19) 


/=“ 4 +"(*) 

To do this, first we reconstruct f{t) from v{t) for a given r*; then we set 
up an optimization procedure to estimate based on the maximization of 

the (negative) /-NS size correlation (a strategy inspired by a similar principle 
was adopted in [^). Fig. panel A, shows an illustrative example of how 
the correlation peaks around the optimal value. As a reference, the dotted 
line marks the value below which 95% of the correlations computed from 
surrogate data fall; surrogate data are obtained by shuffling the values of / 
at the beginning of each NS. 




Figure 7. Slow time-scales inference procedure: test on simulation 
data. Panel A: correlation between low-pass filtered network activity / (see 


Eq. 19) and the size of the immediately subsequent network spike plotted 


against the time-scale r* of the low-pass integrator (continuous line). The 
correlation presents a clear (negative) peak for an ‘optimal’ value = 0.58 
s of the low-pass integrator; such value is interpreted as the effective time-scale 
of the putative slow self-inhibitory mechanism underlying the statistics of 
network events - in this case, short-term synaptic depression (STD); as a 
reference, the dotted line marks the value computed for surrogate data (see 
text). Panel B: for each point in the (wexc, Winh)-plane (see Fig. 
average network activity; the continuous line is the best fit of the theoretical 


expectation for STD’s effective time-scale (Eq. 20); the fitted values for the 


STD parameters tstd and mstd are consistent with the actual values used in 
simulation (tstd = 0.8 s, mstd = 0.2). 


We remark that in this analysis we use both NS and quasi-orbit events 
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(which are both related to the proximity to a Hopf bifurcation). This is 
reasonable since we expect to gain more information about the anti-correlation 
between / and NS size by including both types of large network events. 

Although the procedure successfully recovers a maximum in the correlation, 
the value of (0.58 s) reported in Fig. panel A, is not close to the 

value of tstd (0.8 s). Yet this is expected, since in Eq. 19, t* will in general 


depend on tstd and other parameters of the dynamics, but also on the 
point around which the dynamics is being linearized, more precisely on 
the average activity (u). Specifically, when the fatigue variable follows the 
Tsodyks-Markram model of STD (which of course was actually the case 
in the simulations), linearizing the dynamics of r around a fixed point (r) 
((r) = 1/(1 + wsTD {^) "Tstd)), r behaves as a simple linear integrator with a 
time-constant: 


'r/ptim = t-std (r) = 


"Tstd 


1 + ustd {i^) Tstd 


( 20 ) 


that depends on rgTo, ^std, and (u). 

To test this relationship, we performed the optimization procedure for 
each point of the excitation-inhibition plane. The optimal t* values across 
the excitation-inhibition plane against (u) are plotted in Fi g. panel B 
(dots). The solid line is the best fit of tstd and mstd from Eq. |20|^ which are 
consistent with the actual values used in simulations. 

This result is suggestive of the possibility of estimating from experiments 
the time-scale of an otherwise inaccessible fatigue variable, by modeling it as 
a generic linear integrator, with a “state dependent” time-constant. 

Fig. 1^ shows the outcome of the same inference procedure for two segments 
of experimental recordings. The plot in panel A is qualitatively similar to 
panel A in Fig. although the peak is broader and the maximum correlation 
(in absolute value) is smaller, the t* peak is clearly identified and statistically 
significant (with respect to surrogates, dotted line), thus suggesting a dominant 
time scale for the putative underlying, unobserved fatigue process. However, 
Fig. panel B, clearly shows two significant peaks in the correlation plot; it 
would be natural to interpret this as two fatigue processes, with time scales 
differing by an order of magnitude, simultaneously active in the considered 
recording segment. 

To test the plausibility of this interpretation, we simulated networks with 
simultaneously active STD and spike-frequency adaptation (SFA, see Models 
and Analysis). Fig. shows the results of time scale inference for two cases 
sharing the same time scale for STD (800 ms) and time scale of SFA differing 
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Figure 8. Slow time-scales inference procedure on ex-vivo data. 


Correlation between low-pass filtered network activity / (see Eq. 19) and the 
size of the immediately subsequent network spike plotted against the time- 
scale T* of the low-pass integrator for two experimental datasets (different 
periods - about 40 minutes each - in a long recording). The plot in panel A 
is qualitatively similar to the simulation result shown in panel A of Fig. a 
peak, although broader and of smaller maximum (absolute) value, is clearly 
identified and statistically significant (with respect to surrogate data, dotted 
line). Panel B shows two significant peaks in the correlation plot, a possible 
signature of two concurrently active fatigue processes, with time scales differing 
by roughly an order of magnitude. Panel A: same data as Fig. panel B. 


by a factor of 2 (tsfa = 15 and 30 s respectively). In both cases the negative 
correlation peaks at around t* ~ 500 ms; this peak is plausibly related to 
the characteristic time of STD, consistently with Fig. The peaks at higher 
T*s, found respectively at 12 and 22 s, roughly preserve the ratio of the 
corresponding tsfa values. 

This analysis provides preliminary support to the above interpretation 
of the double peak in Fig. right panel, in terms of two coexisting fatigue 
processes with different time scales. 

We also checked to what extent the avalanche sizes were influenced by 
the immediately preceding amount of available synaptic resources r, and 
we found weak or no correlations; this further supports the interpretation 
offered at the end on the previous section, that avalanches are a genuine 
manifestation of the network excitability which amplifies a wide spectrum of 
small fluctuations. 
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Figure 9. Slow time-scales inference procedure on simulation data 
with STD and spike-frequency adaptation. Correlation between low- 


pass filtered network activity / (see Eq. 19) and the size of the immediately 
snbseqnent network spike plotted against the time-scale t* of the low-pass 
integrator. In this case, the mean-field model includes, besides short-term 
depression (STD), a mechanism mimicking spike-frequency adaptation. Panel 
A: spike-frequency adaptation with characteristic time tsfa = 15 s. Panel B: 
t'sfa = 30 s. In both cases the correlation presents a STD-related peak at 
around r* ~ 500 ms (tstd = 800 ms), consistently with Fig. The peaks at 
higher r*s, found respectively at 11 and 18 s, roughly preserve the ratio of 
the corresponding tsfa values. 


Discussion 


Several works recently advocated a key role of specific network connectivity 
topologies in generating ‘critical’ neural dynamics as manifested in power-law 
distributions of avalanches size and duration (see 


20,^). Also, it has been 


suggested that ‘leader neurons’, or selected coalitions of neurons, play a pivotal 
role in the onset of network events 


(see e.g. 


21 43-45 ). While a role of network 


topology, or heterogeneity in neurons’ excitability, is all to be expected, we 
set out to investigate what repertoire of network events is accessible to 
a network with the simplest, randomly sparse, connectivity, over a wide 
range of excitation-inhibition balance, in the presence of STD as an activity- 
dependent self-inhibition. In the present work we showed that network spikes, 
avalanches and also large fluctuations we termed ‘quasi-orbits’ coexist in such 
networks, with various relative weights and statistical features depending on 
the excitation-inhibition balance, which we explored extensively, including the 
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role of finite-size noise (irregular synchronous regimes in balanced excitatory- 
inhibitory networks has been studied in (^). We remark in passing that 
the occurrence of quasi-orbits is primarily related to the proximity to a Hopf 
bifurcation for the firing rate dynamics; on the other hand, the occurrence of 
NS and, presumably, avalanches, does not necessarily require this condition: 
for instance, NS can occur in the proximity of a saddle-node bifurcation, 
where the low-high-low activity transitions derive from the existence of two 
fixed points, the upper one getting destabilized by the fatigue mechanism 
(see e.g. 46,^); notably, in 12 the authors find that, in a network of leaky 


integrate-and-fire neurons endowed with STD, when a saddle-node separates 
an up- and a down-state, the dynamics develops avalanches during up-state 
intervals only. We also remark that, with respect to the power-law distribution 
of avalanches, it is now widely recognized that while criticality implies power- 
law distributions, the converse is not true, which leaves open the problem of 
understanding what is actually in operation in the neural systems observed 
experimentally (for a general discussion on the issues involved, see j^). In the 
present work, we do not commit ourselves to the issue of whether avalanches 
could be considered as evidence of Self-Organized Criticality. 

In summary, the main contributions of the present work can be listed as 
follows. 

We present a low-dimensional network model, derived from the mean 
field theory for interacting leaky integrate-and-fire neurons with short-term 
depression, in which we include the effect of finite-size (multiplicative) noise. 

At the methodological level we introduce a probabilistic model for events 
detection, and a method for inferring the time-scale(s) of putative fatigue 
mechanisms. At the phenomenological level we recognize the existence of 
quasi-orbits as an additional type of network event, we show the coexistence of 
quasi-orbits, network spikes, and avalanches, and study their different mixing 
depending on the excitability of the network. We also offer a theoretical 
interpretation of the phenomenology, through a bifurcation analysis of the 
mean-field model, and a prediction on the effect of noise in the proximity of 
a Hopf bifurcation. 
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